C.................... MINORA.FOR;1 .......... April 2000
C.... This function is the main sequencing control for the He+ and N+
C.... routines. It is similar to subroutine LOOPS on RSLPSD.FOR. Look 
C.... there for comments.
C.... Cleaned up and commented by P. Richards in April 2000
      SUBROUTINE XION(S,NEQ,NV,TI,ITAU,JBN,SEC,MAXIT,ISOL,IVPLS)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      INTEGER ION,NEQ,NFLAG,LAXIT
      REAL SEC
      DIMENSION N(4,401),TI(3,401),F(20),S(NEQ,1),DCRQ(2,401),V(2)
     > ,DUM(9),NV(4,401)
      COMMON/SLV/DELTA(1806),RHS(1806),WORK(50000)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      REAL SUMION,SUMEXC
      !-- common for total ion production rates EUV+e*. And excitations
      COMMON/BPRTOT/SUMION(3,12,401),SUMEXC(3,12,401)
      COMMON/H/HPLS(401),VHPLS(401)
      COMMON/LPS/SZA(401),EPSN,DC,IMOD,I1,I2,IDNDT,IPMX,IPP,ISKP
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/NEED/HEPLS(401),XPRD(401)
      COMMON/NPRD/COLUM(3,401),OTHPR1(6,401),OTHPR2(6,401),OTHPR3(6,401)
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      COMMON/SAV2/NMSAVE(2,401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
      !.. Switches for wind type and equations solved
      COMMON/CONFAC/IWIND,IVIBN2,IODDN,INPLUS,IHEP,ITEMP,IHPOP
      DATA ICALL/3/,DCR/0.0/,ZLBHE/120.0/,ZLBNP/120.0/
      JB=1
      IPR=1      ! # species: do not change
      IF(IABS(IVPLS).EQ.9)  XMAS=XMASS(3)  !FOR HE+
      IF(IABS(IVPLS).EQ.11) XMAS=XMASS(1)  !FOR N+
      
      ZLBDNP=Z(1)   ! set N+ lower boundary
      IF(ZLBDNP.LT.85) ZLBDNP=85

      !..transferring minor ion densities from storage in XION to N
      DO 41 J=1,JMAX
         !... for N+
         IF(IABS(IVPLS).EQ.11) THEN
            N(1,J)=XIONN(1,J)
            N(3,J)=NV(3,J)-XIONN(1,J)
            IF(N(3,J).LT.0.0) N(3,J)=0.0
            HEPLS(J)=NV(4,J)
            !... N+ production from photoionization
            XPRD(J)=SUMION(3,4,J)+SUMION(3,5,J)+SUMION(3,6,J)
     >         +OTHPR2(3,J)
         ENDIF
         !... for He+
         IF(IABS(IVPLS).EQ.9) THEN
            N(1,J)=NV(4,J)
            N(3,J)=NV(3,J)
            !.. small He+ production added to avoid convergence problems
            HEPRD(J)=OTHPR1(2,J)+1.0E-07
         ENDIF
         !..store minor ion density to calculate dn/dt
         NMSAVE(1,J)=N(1,J)
         !.. store major ion densities and velocities
         OPLS(J)=NV(1,J)
         VOPLS(J)=U(1,J)
         HPLS(J)=NV(2,J)
         VHPLS(J)=U(2,J)
         DCRQ(1,J)=1.0
         DCRQ(2,J)=1.0
         !..  set lower boundary index
         IF(IABS(IVPLS).EQ.9.AND.Z(J).LT.ZLBHE.AND.J.LT.JMAX/2) JB=J
         IF(IABS(IVPLS).EQ.11.AND.Z(J).LT.ZLBNP.AND.J.LT.JMAX/2) JB=J
 41   CONTINUE

      !.. set up key indices for grid
      JBS=JMAX-JB+1
      MIT=JMAX-2*JB+2
      IEQ=IPR*(MIT-2)      !.. number of equations to solve
      MIT1=MIT-1

      !... calculate average values for quantities that don't
      !... change in Newton iteration
      CALL DENAVE(JMAX1,TI,0,0D0,0.0D0)
      CALL AVDEN2(JMAX1,TI,JWR,ZLO,ZHI,IVPLS)
      
      
      !.. Chemical densities for He+ 2020-07-29
      IF(IABS(IVPLS).EQ.9.AND.IHEP.EQ.0) THEN
        DO J=1,JMAX
          CALL MCHEMQ(1,J,JMAX,DT,DUM,N,TI,NMSAVE)
          NV(4,J)=DUM(1)
          DUM(2)=DUM(1)
          IF(DUM(1).LE.0.025*HPLS(J)) NV(4,J)=0.025*HPLS(J)
c          IF(J.EQ.JMAX/2+1) 
c     >       WRITE(40,'(A,I5,2F10.2,1P3E10.2,4X,9E10.2)') ' He+',J,Z(J)
c     >        ,SEC/3600,NV(4,J),DUM(2),0.025*HPLS(J),HPLS(J)
        ENDDO
        RETURN
      ENDIF
      !.. Chemical densities for N+ 2020-07-29
      IF(IABS(IVPLS).EQ.11.AND.INPLUS.EQ.0) THEN
        DO J=1,JMAX
          DUM(1)=0.10*OPLS(J)
          CALL NCHEMQ(1,J,JMAX,DT,DUM,N,TI,NMSAVE)
          XIONN(1,J)=DUM(1)
          DUM(2)=DUM(1)
          IF(DUM(1).GE.0.10*OPLS(J)) XIONN(1,J)=0.10*OPLS(J)
c          IF(J.EQ.JMAX/2+1) 
c     >       WRITE(41,'(A,I5,2F10.2,1P3E10.2,4X,9E10.2)') 
c     >      ' N+',J,Z(J),SEC/3600,XIONN(1,J),DUM(2),0.1*OPLS(J),OPLS(J)
        ENDDO
        RETURN
      ENDIF

      !... Set up max iterations for Newton
      LAXIT=IABS(MAXIT)
      IF(DT.LE.120.AND.LAXIT.LT.29) LAXIT=29

      !***************** Main Loop begins **************************
      DO 220 KTER=1,LAXIT
        !.. set boundary conditions on density
        DO 321 J=1,JMAX
          IF(J.GT.JB.AND.J.LT.JBS) GO TO 321
          IF(IABS(IVPLS).EQ.9) CALL 
     >      MCHEMQ(1,J,JMAX,DT,DUM,N,TI,NMSAVE)
          IF(IABS(IVPLS).EQ.11) CALL NCHEMQ(1,J,JMAX,DT,DUM,N,TI,NMSAVE)
          N(1,J)=DUM(1)
          N(2,J)=DUM(2)
 321    CONTINUE

        !.. call MDFIJ to get unperturbed value to calculate dFij/dn
        DO 145 J=2,MIT1
          KR=IPR*(J-2)
          JC=J+JB-1
          CALL MDFIJ(JC,0,IPR,N,TI,F,JB,JBS,V,IVPLS,XMAS)
          DO 139 IRHS=1,IPR
139       RHS(KR+IRHS)=F(IRHS)
145     CONTINUE
C
        !.. Now set up the Jacobian Matrix dFij/dn
        CALL HMATRX(S,RHS,IEQ,IPR,N,F,JB,JBS,MIT,TI,IVPLS,XMAS)

        !.. invert the jacobian matriX *s* in the inversion routine *bdslv*.
        !.. the increments are stored in array delta in this order
        !.. X(1...n,j),X(1...n,j+1),X(1...n,j+2),....X(1...n,jmaX-1)
        IBW=2*IPR-1
        CALL BDSLV(IEQ,IBW,S,0,RHS,DELTA,WORK,NFLAG)
     
        !.. Check for problems in band solver
        IF(NFLAG.NE.0) THEN
          WRITE(6,*) ' '
          WRITE(6,*) ' *** Problem in band solver ****'
          CALL RUN_ERROR    !.. print ERROR warning in output files
          STOP !3
        ENDIF

        IDIV=0             !... convergence indicator
        DCR=1.0            !... convergence indicator

        !.. test density increments 'dinc' to ensure density > 0 
        !.. (mod steepest descent)
        DO  142 I=1,IPR
          ION=2*IPR-I
          DO 142 J=2,MIT1
            DCRP=1.0
            JC=JB+J-1
            DCRQ(I,JC)=1.0
            DINC=DELTA(IPR*J-ION)
            IF(DINC.LE.0) GO TO 142
            IF(ABS(DINC/N(I,JC)).GT.EPSN) DCRP=DC*ABS(N(I,JC)/DINC)
            IF((DCRP.LT.DCR).AND.(Z(JC).GT.0.0)) DCR=DCRP
            IF(KTER.GT.0) DCRQ(I,JC)=DCRP
 142     CONTINUE
 
        ADCR=DCR
        !.. add iterative increment to the density. 
        !.. convergence when IDIV=0. 
 144    DO 42 I=1,IPR
          ION=2*IPR-I
          DO 42 J=2,MIT1
            JC=JB+J-1
            DINC=DELTA(IPR*J-ION)
            N(I,JC)=N(I,JC)-DINC*ADCR
            IF(ABS(DINC/N(I,JC)).GT.1E-3)  IDIV=IDIV+1
 42     CONTINUE

        !.. if dcr(i) is not close to 1 by the tenth iteration, give up.
        IF(DCR.LE.1.0E-3.AND.KTER.GT.10) GO TO 230

        !.. test to see if convergence has occured.
        IF(IDIV.EQ.0)  GO TO 230
 220  CONTINUE

 230  CONTINUE    !*************** end main loop ***************

      !.. Successful convergence
      !.. Transfer new densities to storage if convergence
      IF(IDIV.EQ.0) THEN
        DO 777 J=1,JMAX
          IF(IABS(IVPLS).EQ.9) NV(4,J)=N(1,J)      ! He+
          IF(IABS(IVPLS).EQ.11) XIONN(1,J)=N(1,J)  ! N+
 777    CONTINUE
      ENDIF

      !.... restore lower boundary to original value after decent time
      IF(SEC-ZLBSEC.GT.3600) ZLBHE=120.0
      IF(SEC-ZLBSEC.GT.3600) ZLBNP=120.0

      IF(IDIV.EQ.0.AND.MAXIT.GT.0) RETURN
      
      !... Unsuccessful - restore all densities to original values
      DO 255 J=1,JMAX
        TI(1,J)=TISAV(1,J)
        TI(2,J)=TISAV(2,J)
        TI(3,J)=TISAV(3,J)
        NV(1,J)=NSAVE(1,J)
        NV(2,J)=NSAVE(2,J)
 255  CONTINUE

      !... raise lower boundary and record the time
      IF(IVPLS.EQ.9.AND.ZLBHE.LT.300) ZLBHE=(ZLBHE+300)/2
      IF(IVPLS.EQ.9)  ZLBSEC=SEC
      IF(IVPLS.EQ.11.AND.ZLBNP.LT.220) ZLBNP=(ZLBNP+220)/2
      IF(IVPLS.EQ.11) ZLBSEC=SEC
      ITER=-50-KTER     !... signals FLIP to reduce time step
      IF(DT.GT.5) RETURN

      !..  non-convergence diagnostics. grid points where problem occured
      IF(IVPLS.EQ.11) THEN
         WRITE(3,116) KTER,DCR
         WRITE(6,116) KTER,DCR
      ELSE
         WRITE(3,117) KTER,DCR
         WRITE(6,117) KTER,DCR
      ENDIF
 116    FORMAT(/'  NON-CONVERGENCE IN N+ SOLUTION, ITER='
     > ,I3,'  DCR1=',1P,E10.1)
 117    FORMAT(/'  NON-CONVERGENCE IN He+ SOLUTION, ITER='
     > ,I3,'  DCR1=',1P,E10.1)


      !... Write out problem points on the field line
      DO 535 J=JB+1,JBS-1
        IF(DCRQ(1,J)+DCRQ(2,J).LT.2) THEN
          IF(J.LE.JMAX/2) THEN
            WRITE(6,537) J,INT(Z(J)),DCRQ(1,J),DCRQ(2,J)
          ELSE
            WRITE(6,538) J,INT(Z(J)),DCRQ(1,J),DCRQ(2,J)
          ENDIF
        ENDIF
 535  CONTINUE
 537  FORMAT(2X,'Northern hemisphere:- problem at point=',I5,' ALT=',I8
     > ,1P,4E8.1)
 538  FORMAT(2X,'Southern hemisphere:- problem at point=',I5,' ALT=',I8
     > ,1P,4E8.1)

      !.... Return ITER=-1 to stop execution
      ITER=-1
      RETURN
      END
C::::::::::::::::::::::::::: MCHEMQ :::::::::::::::::::::::::::::::::::::
C... Gets the chemical equilibrium density for He+
      SUBROUTINE MCHEMQ(ISW,J,JMAX,DT,DUM,N,TI,NMSAVE)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      DIMENSION RTS(199),N(4,401),TI(3,401),NMSAVE(2,401),DUM(9)
      CALL RATS(J,TI(3,J),TI(1,J),TN(J),RTS)
        !.. 100000.0 /cc prevents divide problems when called at high altitudes
        HELOSS=(RTS(44)+RTS(45))*(N2N(J)+100000.0)+
     >    (RTS(75)+RTS(76))*O2N(J)
        DUM(1)=HEPRD(J)/HELOSS   !... production/loss
        DUM(2)=0.0
      RETURN
      END
C::::::::::::::::::::::::::::: NCHEMQ :::::::::::::::::::::::::::::::::::::
C... Gets the chemical equilibrium density for N+
      SUBROUTINE NCHEMQ(ISW,J,JMAX,DT,DUM,N,TI,NMSAVE)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/NEED/HEPLS(401),XPRD(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/MIB/OP2D(401),OP2P(401),N2D(401),N4S(401),NNO(401)
      DIMENSION RTS(199),N(4,401),TI(3,401),NMSAVE(2,401),DUM(9)
        CALL RATS(J,TI(3,J),TI(1,J),TN(J),RTS)
        !... Evaluate loss rates
        QNP1=RTS(30)*O2N(J)
        QNP2=RTS(25)*O2N(J)
        QNP3=RTS(59)*O2N(J)
        QNP4=RTS(31)*ON(J)
        QNP5=RTS(22)*O2N(J)
        QNP6=RTS(65)*O2N(J)
        QNP7=RTS(66)*O2N(J)
        !.. 1.00 /cc prevents divide problems when called at high altitudes
        TLNPLUS=RTS(31)*ON(J)+
     >   (RTS(22)+RTS(25)+RTS(30)+RTS(59)+RTS(65)+RTS(66))*(O2N(J)+1.00)
        !... evaluate production rates
        PNP1=RTS(45)*HEPLS(J)*N2N(J)
        PNP2=XPRD(J)
        PNP3=RTS(29)*N2D(J)*OPLS(J)
        PNP4=RTS(22)*OP2P(J)*N2N(J)*0.000  !.. does not happen

        !.. chemical equilibrium density for N+
        DUM(1)=(PNP1+PNP2+PNP3+PNP4)/TLNPLUS
        DUM(2)=0.0
        !...IF(ISW.EQ.1) WRITE(35,*)PNP1,PNP2,PNP3,PNP4,QNP1,QNP2
        !...> ,QNP3,QNP4
      RETURN
      END
C::::::::::::::::::::::::::::: MDFIJ ::::::::::::::::::::::::::::::::::::
C.... This routine sets up the He+ and H+ continuity equations to be 
C.... minimized by the Newton solver. See program DFIJ on RSDENA.FOR 
C.... for more details

      SUBROUTINE MDFIJ(J,JSJ,IPR,N,TI,F,JB,JBS,V,IVPLS,XMAS)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      REAL MASS
      INTEGER ION
      DIMENSION VL(2),VU(2),FLU(2),FLL(2),V(2),ANM(2),PM(2),Q(2),L(2)
     > ,TINCR(2),N(4,401),TI(3,401),F(20),DIVEM(2)
      COMMON/MSS2/MASS(4),A(4)
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/SAV2/NMSAVE(2,401)
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
      DATA FLU/0.0,0.0/

      !... ID is the unit number for writing diagnostics in
      !... both hemispheres
      ID=3  
      IF(J.GT.JMAX/2) ID=9
      IF(JSJ.NE.9) ID=0

      !.. MINVEL to get the velocities(fluxes) at the lower 1/2 pt
      FLL(1)=FLU(1)
      FLL(2)=FLU(2)
      CALL MINVEL(J-1,VL,FLL,N,JSJ,0,IPR,XMAS,IVPLS) ! f90 TI not referenced

      !..  CALL CHEMO to evaluate the source (Q) and sink(L) terms 
      IF(IABS(IVPLS).EQ.9) CALL CHEMO9(IPR,J,Q,N,NMSAVE,TI,1,L,DIVEM)
      IF(IABS(IVPLS).EQ.11) CALL CHEM11(IPR,J,Q,N,NMSAVE,TI,1,L,DIVEM)

      !... call dave for dn/dt. PM= dn(t+dt)/dt, ANM= dn(t)/dt
      CALL DAVEM(ITER,IPR,J,ANM,PM,N,NMSAVE,BM)

      !.. MINVEL to get the velocities(fluxes) at the upper 1/2 pt
360   CALL MINVEL(J,VU,FLU,N,JSJ,ID,IPR,XMAS,IVPLS)

      !.. magnetic field strength at midpoints
      BU=(BM(J)+BM(J+1))*0.5
      BL=(BM(J)+BM(J-1))*0.5

      !.. Set up F=prod-loss-dn/dt-div(flux)
      DO 380 I=1,ION
      TINCR(I)=(PM(I)-ANM(I))/DT     !.. dn/dt
      FGR=(FLU(I)/BU-FLL(I)/BL)      !.. div(flux)
C
      F(I)=Q(I)-L(I)-TINCR(I)-FGR-DIVEM(I)

      !.... average velocity for He+ and N+ ........
      IF(JSJ.EQ.0.AND.IABS(IVPLS).EQ.9) UE(I,J)=.5*(VL(I)+VU(I))
      IF(JSJ.EQ.0.AND.IABS(IVPLS).EQ.11) XIONV(1,J)=.5*(VL(I)+VU(I))
      VL(I)=VU(I)

      !.. printing of factors in the continuity equations ...........
      !.. since all quantities have been divided by bm and integrated,
      !.. multiply by bd to get actual magnitudes
      IF(ID.EQ.0) GO TO 380
        BD=2.0*BM(J)/(SL(J+1)-SL(J-1))
        WRITE(ID,666) J,BD,FLU(I),FLL(I),FGR,Q(I),L(I),F(I),UE(I,J)
     >   ,TINCR(I),XIONV(1,J),DIVEM(I)
  666    FORMAT(2X,'FIJ2',I4,1P,22E11.3)
         IF(I.EQ.2) WRITE(ID,*) '   '
  380 CONTINUE
      RETURN
      END
C::::::::::::::::::::::::::::: MINVEL ::::::::::::::::::::::::::::::::::
C... This subroutine returns the ion fluxes for the continuity equation.
      SUBROUTINE MINVEL(J,V,FLUX,N,JSJ,ID,IPR,XMAS,IVPLS)  ! f90 added
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      REAL MASS
      DIMENSION PP(4),QSIGN(2),V(2),FLUX(2),GAMMA(2),B(3)
      DIMENSION N(4,401),HX(2)
      COMMON/OTD/OTDC(2)
      COMMON/H/HPLS(401),VHPLS(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      COMMON/MSS2/MASS(4),A(4)
      COMMON/NEED/HEPLS(401),XPRD(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/SAV2/NMSAVE(2,401)
      COMMON/GRAD/GRADTE(401),GRADTI(401)
      DATA QSIGN/-1.,1./

      !... evaluate total ion density
      NEU=OPLS(J+1)+HPLS(J+1)
      NEL=OPLS(J)+HPLS(J)
      NE=OP(J)+DSQRT(HPLS(J)*HPLS(J+1))
      IF(IABS(IVPLS).GT.9) THEN
         NEU=NEU+HEPLS(J+1)+N(1,J+1)+N(3,J+1)
         NEL=NEL+HEPLS(J)+N(1,J)+N(3,J)
         NE=NE+DSQRT(HEPLS(J)*HEPLS(J+1))+DSQRT(N(1,J+1)*N(1,J))
         NE=NE+DSQRT(N(3,J)*N(3,J+1))
      ENDIF
      !.. interpolate minor ion to midpoint
      PP(1)=DSQRT(N(1,J+1)*N(1,J))
      PP(2)=0.0
      !.. CALCULATE ELECTRON DENSITY AT UPPER, LOWER AND MID POINT
      IF (IABS(IVPLS).LE.9) THEN
        NEU=NEU+N(1,J+1)+N(3,J+1)
        NEL=NEL+N(1,J)+N(3,J)
        NE=NE+PP(1)+DSQRT(N(3,J)*N(3,J+1))
      ENDIF

      !.. CALL DCOEFF to obtain d(i),hX,beta1,beta2,betaX, @ b(i) used
      !.. to calc v(i)  all from st. maurice and schunk
      CALL DCOEFF(J,PP,D,B,HX,GAMMA,JSJ,ID,SUMN,XMAS)    ! f90 added

      !-----calculate altitude derivatives----------
      !.. note that TEJ(J)=0.5(TI(3,J+1)+TI(3,J))/DS(J)/TIJ(J) see AVDEN
      GRADNE=TEJ(J)*(NEU-NEL)/NE

      DLNI=(N(1,J+1)-N(1,J))/DS(J)/PP(1)              !.. dn/ds
      GRA=XMAS*GR(1,J)                                !.. gravity
      QION=QSIGN(2)*GRADTE(J)

      !.... neutral wind term SUMN=sum of collision terms for ions on neutrals
      Q2=UNJ(J)*SUMN
      !.. Q3 is the sum of FRICTION terms.  HX(1)=X - H+, HX(2)=X - O+
      Q3=0.5*(VOPLS(J+1)+VOPLS(J))*HX(2)
      Q3=Q2+Q3+0.5*(VHPLS(J+1)+VHPLS(J))*HX(1)
      Q4=GRADTI(J)*(1.0D00-(B(1)+B(2)+B(3)))

      !.. Velocity from the momentum equation.
       V(1)=Q3-D*(DLNI-GRA+Q4+GRADNE+QION)

      !.......  print routine for momtm eqn terms ..........
      IF(ID.NE.0) WRITE(ID,605) J,Z(J),DLNI,GRADTI(J),GRADNE
     > ,GRADTE(J),GRA,QION,Q4,Q2,Q3
 605   FORMAT(1X,'MINVEL',I4,F9.0,1P,22E11.3)

       !.. ion flux
       FLUX(1)=V(1)*PP(1)
       !... print velocities ...........
      IF(JSJ.EQ.-4) WRITE(3,606) J,Z(J),V(1),V(2)
 606  FORMAT(1X,'J=',I4,3X,'ALT=',F9.0,2X,'V1=',1P,E13.6,2X,'V2=',E13.6)
      RETURN
      END
C::::::::::::::::::::::::::::: DCOEFF ::::::::::::::::::::::::::::::::::::
C... This routine calculates the diffusion velocities from the minor ion
C... equation 34 as given by St. Maurice and Schunk Planet and Spac. 
C... Sci. 1976.  Units are cgs
      SUBROUTINE DCOEFF(J,PP,D,B,HX,GAMMA,JSJ,ID,SUMN,XMAS)
      IMPLICIT DOUBLE PRECISION (A-H,M-Z)
      REAL*8 NU(3,3),NUPI(3,2),N(3),NUP(3),MU(3,3),KB
      INTEGER MM
      DIMENSION D1(3,3),D4(3,3),A(3),AM(3),PP(4)
      DIMENSION B(3),DELT(2),HX(2),R(2,2),Y(3,2),GAMMA(2)
      COMMON/H/HPLS(401),VHPLS(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      DATA A/1.0D00,16.0D00,0.0D00/,KB/1.38062D-16/
      DATA AMASS/1.67D-24/,RE/6.371D08/

      !.. set up constants: A= amu of mass A1=hydrogen,A2=oxygen,
      !.. C A3=helium. amass=1.67e-24 grams
      AM(1)=A(1)*AMASS
      AM(2)=A(2)*AMASS
      AM(3)=XMAS
      A(3)=XMAS/AMASS
      !.. load densities H+,O+,and minor ions.
      N(1)=DSQRT(HPLS(J+1)*HPLS(J))
      N(2)=OP(J)
      N(3)=PP(1)
      !..  Mass dependent constants
      DO I=1,3
        AI2=A(I)*A(I)
          DO K=1,3
            IF (K.EQ.3.OR.K.EQ.I) GO TO 11
            AK2=A(K)*A(K)
            AIK2=(A(I)+A(K))*(A(I)+A(K))
            D1(I,K)=(3.0*AI2+.1*A(I)*A(K)-.2*AK2)/AIK2
            D4(I,K)=(1.2*AK2-1.5*A(I)*A(K))/AIK2
11          MU(I,K)=A(I)*A(K)*AMASS/(A(I)+A(K))
          END DO
      END DO

      !.. Evaluate collision terms. Neutral terms from AVDEN2
      SUMN=NUX(1,J)   !.. sum of neutral collision terms

      DO I=1,3
        DO K=1,2
          AST=A(I)*A(K)/(A(I)+A(K))
          NU(I,K)=1.27D00*DSQRT(AST)*N(K)/(A(I)*TIJ(J)**1.5)
        END DO
      END DO

      NUP(1)=NU(1,1)+1.25*NU(1,2)*(D1(1,2)+1.5*MU(1,2)/AM(1))
      NUP(2)=NU(2,2)+1.25*NU(2,1)*(D1(2,1)+1.5*MU(2,1)/AM(2))
      NUP(3)=1.25*NU(3,1)*(D1(3,1)+1.5*MU(3,1)/AM(3))
      NUP(3)=NUP(3)+1.25*NU(3,2)*(D1(3,2)+1.5*MU(3,2)/AM(3))
      NUPI(1,2)=1.25*NU(1,2)*(D4(1,2)+1.5*MU(1,2)/AM(1))
      NUPI(2,1)=1.25*NU(2,1)*(D4(2,1)+1.5*MU(2,1)/AM(2))
      NUPI(3,1)=1.25*NU(3,1)*(D4(3,1)+1.5*MU(3,1)/AM(3))
      NUPI(3,2)=1.25*NU(3,2)*(D4(3,2)+1.5*MU(3,2)/AM(3))
      EPS=1.0/(1.0-(NUPI(1,2)/NUP(1))*(NUPI(2,1)/NUP(2)))

      DO I=1,2
        DO K=1,2
          IF (K.EQ.I) THEN
          MM=3-K
          Y(I,K)=1.0-NUPI(3,I)/NUP(3)-NUPI(MM,I)/NUP(MM)*
     >          NUPI(3,MM)/NUP(3)
          ELSE
          Y(I,K)=NUPI(K,I)/NUP(K)-NUPI(3,I)/NUP(3)-
     >          NUPI(K,I)/NUP(K)*NUPI(3,K)/NUP(3)
          ENDIF
        END DO
      END DO

      Y(3,1)=-1.0/EPS
      Y(3,2)=Y(3,1)
      X1=MU(3,1)/AM(3)*NU(3,1)/NUP(3)+MU(3,2)/AM(3)*NU(3,2)/NUP(3)
      DELT(1)=9.0*(MU(3,1)/AM(3))*X1/8.0
      DELT(2)=9.0*(MU(3,2)/AM(3))*X1/8.0
      DTOP=KB*TIJ(J)/AM(3)
      NUSUM=NU(3,1)*(1.0-DELT(1))+NU(3,2)*(1.0-DELT(2))+SUMN

      !.. divide diffusion and neutral collision terms by total of
      !.. of collision terms nusum. redefine neutral collision term sumn
      D=DTOP/NUSUM
      SUMN=SUMN/NUSUM
      BTOPC=15.0*AM(3)*EPS/8.0
      B(1)=BTOPC/AM(1)*(MU(3,1)/AM(1)*NU(3,1)/NUP(1)*Y(1,1)
     >       +MU(3,2)/AM(1)*NU(3,2)/NUP(1)*Y(1,2))
      B(2)=BTOPC/AM(2)*(MU(3,1)/AM(2)*NU(3,1)/NUP(2)*Y(2,1)
     >       +MU(3,2)/AM(2)*NU(3,2)/NUP(2)*Y(2,2))
      B(3)=BTOPC/AM(3)*(MU(3,1)/AM(3)*NU(3,1)/NUP(3)*Y(3,1)
     >       +MU(3,2)/AM(3)*NU(3,2)/NUP(3)*Y(3,2))

      !...  UNEUT=SUMN*VNEUT*1.0D02/NUSUM
      DO I=1,2
        DO K=1,2
          IF (K.NE.I) THEN
            R(I,K)=MU(I,K)/AM(I)*NU(I,K)/AM(I)*EPS/NUP(I)
            X1=MU(3,I)*NU(3,I)*Y(I,I)+MU(3,K)*NU(3,K)*Y(I,K)
            R(I,K)=9.0*R(I,K)*X1/8.0
          ENDIF
        ENDDO
      END DO

      !..HX=velocity coefficients (friction)
        HX(1)=(NU(3,1)*(1.0-DELT(1))-R(1,2)+R(2,1))/NUSUM
        HX(2)=(NU(3,2)*(1.0-DELT(2))-R(2,1)+R(1,2))/NUSUM
      RETURN
      END
C:::::::::::::::::::::::::::::::: CHEMO9 :::::::::::::::::::::::::::::::
C.... this program determines the interpolated production and loss
C.... processes for He+. It CALLs rates to get the rate constants and TERLIN
C.... to do the interpolation. Updated July 1998 to include ExB drift
C.... See RSDENA_ExB.F for details. Bailey et al. PSS 1983,page 389 (DIVEM)
      SUBROUTINE CHEMO9(ION,JI,SOURCE,N,NMSAVE,TI,ISW,SINK,DIVEM)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION N(4,401),TI(3,401),NMSAVE(2,401),Q(2,3),L(2,3)
     > ,SOURCE(2),SINK(2),RTS(199),DIVEM(2),DV(2,3)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/ALT/Z(401),DT,DH,FHT,EPS,TF,ITF,JMAX,JMAX1,ITER,JION
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      COMMON/OHPL/OLOSS(401),HLOSS(401),HPROD(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      common/ExB/eq_vperp,div_vem(401)

      !-- Q and L are chemical source and sink, DV is for ExB drift
      DO 300 K=1,3
        J=K+JI-2
        CALL RATS(J,TI(3,J),TI(1,J),TN(J),RTS)
        HELOSS=(RTS(44)+RTS(45))*N2N(J)+(RTS(75)+RTS(76))*O2N(J)
        HPRODN=HPROD(J)*OPLS(J)
        Q(1,K)=HEPRD(J)*RBM(J)
        Q(2,K)=HPRODN*RBM(J)
        L(1,K)=HELOSS*N(1,J)*RBM(J)
        L(2,K)=HLOSS(J)*N(2,J)*RBM(J)
        !.. Silvy ExB drift term is added. dv(1,k) is for He+/N+ and dv(2,k) for H+
        DV(1,K)=div_vem(J)*(0.5*(N(1,J)+NMSAVE(1,J)))*RBM(J)
        DV(2,K)=div_vem(J)*(0.5*(N(2,J)+NMSAVE(2,J)))*RBM(J)
 300  CONTINUE

      !.... terd is called to interpolate
      DO 400 I=1,ION
        CALL TERLIN(JI,Q(I,1),Q(I,2),Q(I,3),DS(JI-1),DS(JI),QL,QM,QU)
        CALL TERLIN(JI,L(I,1),L(I,2),L(I,3),DS(JI-1),DS(JI),LL,LM,LU)
        CALL TERLIN(JI,DV(I,1),DV(I,2),DV(I,3),DS(JI-1),DS(JI),DL,DM,DU)
        SOURCE(I)=QM
        SINK(I)=LM
        DIVEM(I)=DM*eq_vperp
 400  CONTINUE
      RETURN
      END
C:::::::::::::::::::::::::: CHEM11 ::::::::::::::::::::::::::::::::::::::
C.... this program determines the interpolated production and loss
C.... processes for N+. It CALLs rates to get the rate constants and TERLIN
C.... to do the interpolation.Updated July 1998 to include ExB drift
C.... See RSDENA_ExB.F for details. Bailey et al. PSS 1983,page 389 (DIVEM)
      SUBROUTINE CHEM11(ION,JI,SOURCE,N,NMSAVE,TI,ISW,SINK,DIVEM)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION N(4,401),TI(3,401),NMSAVE(2,401),Q(2,3),L(2,3)
     > ,SOURCE(2),SINK(2),RTS(199),DIVEM(2),DV(2,3)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      COMMON/NEED/HEPLS(401),XPRD(401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/MIB/OP2D(401),OP2P(401),N2D(401),N4S(401),NNO(401)
      common/ExB/eq_vperp,div_vem(401)

      DO 300 K=1,3
        J=K+JI-2
        CALL RATS(J,TI(3,J),TI(1,J),TN(J),RTS)
        PNP1=RTS(45)*HEPLS(J)*N2N(J)
        PNP2=XPRD(J)
        PNP3=RTS(29)*N2D(J)*OPLS(J)
        PNP4=RTS(22)*OP2P(J)*N2N(J)*0.000   !. does not happen
        LNP1=RTS(30)*O2N(J)
        LNP2=RTS(25)*O2N(J)
        LNP3=RTS(59)*O2N(J)
        LNP4=RTS(31)*ON(J)
        LNP5=RTS(22)*O2N(J)
        LNP6=RTS(65)*O2N(J)
        LNP7=RTS(66)*O2N(J)

        Q(1,K)=(PNP1+PNP2+PNP3+PNP4)*RBM(J)
        Q(2,K)=0.0
        L(1,K)=(LNP1+LNP2+LNP3+LNP4+LNP5+LNP6+LNP7)*N(1,J)*RBM(J)
        L(2,K)=0.0
        !.. Silvy ExB drift term is added. dv(1,k) is for He+/N+ and dv(2,k) for H+
        DV(1,K)=div_vem(J)*N(1,J)*RBM(J)
        DV(2,K)=div_vem(J)*N(2,J)*RBM(J)
 300  CONTINUE
      !... terd is called to interpolate
      DO 400 I=1,ION
        CALL TERLIN(JI,Q(I,1),Q(I,2),Q(I,3),DS(JI-1),DS(JI),QL,QM,QU)
        CALL TERLIN(JI,L(I,1),L(I,2),L(I,3),DS(JI-1),DS(JI),LL,LM,LU)
        CALL TERLIN(JI,DV(I,1),DV(I,2),DV(I,3),DS(JI-1),DS(JI),DL,DM,DU)
        SOURCE(I)=QM
        SINK(I)=LM
        DIVEM(I)=DM*eq_vperp
 400  CONTINUE
      RETURN
      END
C::::::::::::::::::::::::::::: AVDEN2 ::::::::::::::::::::::::::::::::::::
C.....................................................................
C   this program evaluates the interpolated densities at the midpoints
C   it also evaluates the ion-neutral collision frequencies nuX(i,j)
C   and o+ and h+ production and loss rates
C........................................................................
      SUBROUTINE AVDEN2(JMAX1,TI,JWR,ZLO,ZHI,IVPLS)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      REAL OA,HA,N2A,O2A,TNJ,TR,HEA
      DIMENSION TI(3,401)
      COMMON/A/OPLS(401),VOPLS(401),HE(401),UE(2,401),HEPRD(401),OP(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      COMMON/GRAD/GRADTE(401),GRADTI(401)
      COMMON/OHPL/OLOSS(401),HLOSS(401),HPROD(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX2,ITER,ION

      !.. ion-neutral coll freqs taken from schunk and nagy rev. geophys. v18
      !.. p813, 1980. first,  resonant charge eXchange from table 5 p823
      DO 10 J=1,JMAX1
        OP(J)=DSQRT(OPLS(J+1)*OPLS(J))
        OA=SQRT(ON(J)*ON(J+1))
        HA=SQRT(HN(J)*HN(J+1))
        N2A=SQRT(N2N(J)*N2N(J+1))
        O2A=SQRT(O2N(J)*O2N(J+1))
        HEA=SQRT(HE(J)*HE(J+1))
        TNJ=0.5*(TN(J)+TN(J+1))
        TR=(TIJ(J)+TNJ)*0.5
        IF(IABS(IVPLS).EQ.9) CALL ADEN9(J,OA,HA,N2A,O2A,HEA,TIJJ
     >  ,TR,TNJ,NUX)
        IF(IABS(IVPLS).EQ.11) CALL ADEN11(J,OA,HA,N2A,O2A,HEA
     >  ,TIJJ,TR,TNJ,NUX)
10    CONTINUE
      RETURN
      END
C:::::::::::::::::: ADEN9::::::::::::::::::::::::::::::
C.........For HE+
      SUBROUTINE ADEN9(J,OA,HA,N2A,O2A,HEA,TIJJ,TR,TNJ,NUX)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      REAL OA,HA,N2A,O2A,TNJ,TR,SQT,HEA
      DIMENSION NUX(2,401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX2,ITER,ION
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      SQT=SQRT(TR)
      RHEPHE=8.73E-11*HEA*SQT*(1-0.093*ALOG10(TR))**2

      !.. non-resonant ion neutral interactions table 6
      CHEPN=(4.71*HA+10.1*OA+16*N2A+15.3*O2A)*1.0E-10
      NUX(1,J)=RHEPHE+CHEPN
      NUX(2,J)=0.0

      !... print diagnostics
      IF(JWR.NE.2.OR.Z(J).LT.ZLO.OR.Z(J).GT.ZHI) GO TO 10
        ID=3
        IF(J.GT.JMAX/2) ID=9
        WRITE(ID,99) J,Z(J),OA,HA,N2A,O2A,HEA,TIJJ,TNJ,TR,RHEPHE
        WRITE(ID,99) J,Z(J),CHEPN,NUX(1,J),NUX(2,J)
 10   CONTINUE
 99   FORMAT(1X,'AVDEN9',I4,F7.0,1P,22E10.3)
      RETURN
      END
C:::::::::::::::::: ADEN11::::::::::::::::::::::::::::::
C........For N+
      SUBROUTINE ADEN11(J,OA,HA,N2A,O2A,HEA,TIJJ,TR,TNJ,NUX)
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      REAL OA,HA,N2A,O2A,TNJ,TR,SQT,HEA
      DIMENSION NUX(2,401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX2,ITER,ION
      COMMON/PALT/ZLO,ZHI,JWR,JTI
      SQT=SQRT(TR)

      !..  resonant ion neutral collision. No neutral N at present(8/17/89)
      NNA=0.0
      CRIN=3.83D-11*NNA*SQT*(1.0-0.063*ALOG10(TR))**2
      CNRIN=(1.45*HA+1.49*HEA+4.42*OA+7.47*N2A+7.25*O2A)
      CNRIN=CNRIN*1.0D-10
      NUX(1,J)=CRIN+CNRIN
      NUX(2,J)=0.0

      !... print diagnostics
      IF(JWR.NE.2.OR.Z(J).LT.ZLO.OR.Z(J).GT.ZHI) GO TO 10
        ID=3
        IF(J.GT.JMAX/2) ID=9
        WRITE(ID,99) J,Z(J),OA,HA,N2A,O2A,HEA,TIJJ,TNJ,TR,CRIN
        WRITE(ID,99) J,Z(J),CNRIN,NUX(1,J),NUX(2,J)
 10   CONTINUE
 99   FORMAT(1X,'AVDEN11',I4,F7.0,1P,22E10.3)
      RETURN
      END
C::::::::::::::::::::::::::::: HMATRX ::::::::::::::::::::::::::::::
C......  HMATRX EVALUATES DELF/DELT USING STEPHANSONS METHOD.
C......  S(L,M) = DEL(FIJ) / DEL(N)
C......  WHERE N IS THE DENSITY OF ION IV AT ALTITUDE JV. L AND M
C......  ARE COMPUTED FROM JF...IV. THE ARRAY RHS CONTAINS VALUES OF
C...... FIJ SAVED FROM PREVIOUS COMPUTATION
      SUBROUTINE HMATRX(S,RHS,INEQ,IPR,N,F,JB,JBS,MIT,TI,IFLAG,XMA)
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      DIMENSION RHS(1806),S(INEQ,8),N(4,401),TI(3,401),F(20),V(2)
      IBW=4*IPR-1    !.... band width

      DO 190 JZS=1,INEQ
      DO 180 KZS=1,IBW
 180  S(JZS,KZS)=0.0
 190  CONTINUE

      DO 80 JF=2,MIT-1
        J1=MAX0(2,JF-1)
        J2=MIN0(JF+1,MIT-1)
        DO 90 JV=J1,J2
          DO 130 IV=1,IPR
            L=IPR*(JF-2)
            IF(JF.LE.3)M=IPR*(JV-2)+IV
            IF(JF.GT.3)M=IPR*(JV-JF)+IV+2*IPR
            JVC=JB+JV-1
            JFC=JB+JF-1
            H=1.E-4*N(IV,JVC)
            N(IV,JVC)=N(IV,JVC)+H
            CALL MDFIJ(JFC,1,IPR,N,TI,F,JB,JBS,V,IFLAG,XMA)
            N(IV,JVC)=N(IV,JVC)-H
            RH=1/H
            DO 75 IS=1,IPR
              IF(JF.LE.3) S(L+IS,M)=(F(IS)-RHS(L+IS))*RH
              IF(JF.GT.3) S(L+IS,M-IS)=(F(IS)-RHS(L+IS))*RH
 75         CONTINUE
 130      CONTINUE
  90    CONTINUE
  80  CONTINUE
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE DENAVE(JMAX1,TI,JSJ,ZLO,ZHI)
C.....................................................................
C   this program evaluates the interpolated densities at the midpoints
C   it also evaluates the ion-neutral collision frequencies nux(i,j)
C   and o+ and h+ production and loss rates
C........................................................................
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      DIMENSION TI(3,401),RTS(199)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      COMMON/GRAD/GRADTE(401),GRADTI(401)
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      COMMON/SAV2/NMSAVE(2,401)
      COMMON/OHPL/OLOSS(401),HLOSS(401),HPROD(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX2,ITER,ION
      COMMON/CFACT/DHDU_FAC,COLFAC,CFAC(9),ISCH(9)
      COMMON/RK/VC(401)
      DATA BK/1.3807E-16/
      DO 10 J=1,JMAX1
      RBM(J)=1.0/BM(J)
      DS(J)=SL(J+1)-SL(J)
      TIJ(J)=0.5*(TI(1,J)+TI(1,J+1))
      TEJ(J)=0.5*(TI(3,J)+TI(3,J+1))/DS(J)/TIJ(J)
      GR(1,J)=0.5*(GR(2,J)+GR(2,J+1))/BK/TIJ(J)
      GRADTI(J)=(TI(1,J+1)-TI(1,J))/DS(J)/TIJ(J)
      GRADTE(J)=(TI(3,J+1)-TI(3,J))/DS(J)/TIJ(J)
      UNJ(J)=0.5*(UN(J+1)+UN(J))

      !... ion-neutral coll freqs taken from schunk and nagy rev. geophys. v18
      !... p813, 1980. first,  resonant charge exchange from table 5 p823
      OA=SQRT(ON(J)*ON(J+1))
      HA=SQRT(HN(J)*HN(J+1))
      N2A=SQRT(N2N(J)*N2N(J+1))
      O2A=SQRT(O2N(J)*O2N(J+1))
      TNJ=0.5*(TN(J)+TN(J+1))
      TR=(TIJ(J)+TNJ)*0.5
      SQT=SQRT(TR)
      RHPH=2.65E-10*HA*SQT*(1-0.083*DLOG10(TR))**2

      !.. COLFAC is the scaling factor for O+ - O collision frequency  
      !.. should be 1.7 according to Burnside 1987 (in press)
      ROPO= COLFAC * 3.67E-11*OA*SQT*(1-0.064*DLOG10(TR))**2 
      RHPO=6.61E-11*OA*SQRT(TIJ(J))*(1-0.047*DLOG10(TIJ(J)))**2
      ROPH=4.63E-12*HA*SQRT(TNJ+TIJ(J)/16.0)

      !.. non-resonant ion neutral interactions table 6
      CHPN=(33.6*N2A+32*O2A)*1.0E-10
      COPN=(6.28*N2A+6.64*O2A)*1.0E-10
      NUX(1,J)=ROPO+ROPH+COPN
      NUX(2,J)=RHPH+RHPO+CHPN

      !..... O+ and H+ production and loss rates
      CALL RATS(J,TI(3,J),TI(1,J),TN(J),RTS)
      !... multiply O+ + N2 reaction rate by N2* factor
      RTS(3)=RTS(3)*VC(J)
      OLOSS(J)=RTS(2)*HN(J)+RTS(3)*N2N(J)+RTS(4)*O2N(J)
      HLOSS(J)=RTS(1)*ON(J)
      HPROD(J)=RTS(2)*HN(J)

      !... printing results
      IF(JSJ.NE.4.OR.Z(J).LT.ZLO.OR.Z(J).GT.ZHI) GO TO 10
        IUN=34
        !... IF(J.GT.JMAX/2) IUN=9
        WRITE(IUN,99) J,Z(J),OA,HA,N2A,O2A,TIJ(J),TNJ,TR,RHPH,ROPO,RHPO
        WRITE(IUN,99) J,Z(J),ROPH,CHPN,COPN,NUX(1,J),NUX(2,J),OLOSS(J)
     >   ,HLOSS(J),HPROD(J)
 10   CONTINUE
 99   FORMAT(1X,'AVDEN',I4,F7.0,1P,22E10.3)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C.... calculates ante and post values of half interval, and average
C.... of (ion density)/(magnetic field)
      SUBROUTINE DAVEM(ITER,ION,J,ANM,PM,N,NMSAVE,BM)
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,JITER,JION
      DIMENSION ANM(2),PM(2),N(4,401),NMSAVE(2,401),BM(401)

      DO 100 I=1,ION
        !...     ante values
        BB=NMSAVE(I,J-1)*RBM(J-1)
        C=NMSAVE(I,J)*RBM(J)
        A=NMSAVE(I,J+1)*RBM(J+1)
        IF(A*BB*C.LT.0) WRITE(6,115) I,J,ITER,A,BB,C
 115    FORMAT('-VE ANTE-> I,J,ITER,A,B,C',3I6,1P,9E11.2)
        JEQ = JMAX/2+1
        IF(J.LT.JEQ) CALL TERMID(J,BB,C,A,DS(J-1),DS(J),QL,QM,QU)
        IF(J.GE.JEQ) CALL TERMID(J,A,C,BB,DS(J),DS(J-1),QL,QM,QU)
        ANM(I)=QM
        !... post values
        BB=N(I,J-1)*RBM(J-1)
        C=N(I,J)*RBM(J)
        A=N(I,J+1)*RBM(J+1) 
        IF(A*BB*C.LT.0) WRITE(6,116) I,J,ITER,A,BB,C
 116    FORMAT('-VE POST-> I,J,ITER,A,B,C',3I6,1P,9E11.2)
        IF(J.LT.JEQ) CALL TERMID(J,BB,C,A,DS(J-1),DS(J),QL,QM,QU)
        IF(J.GE.JEQ) CALL TERMID(J,A,C,BB,DS(J),DS(J-1),QL,QM,QU)
        PM(I)=QM
 100  CONTINUE
      RETURN
      END
C::::::::::::::::::::::::::: MINBLK :::::::::::::::::::
      BLOCK DATA MINBLK
      IMPLICIT DOUBLE PRECISION(A-H,L,N-Z)
      REAL MASS
      COMMON/MINOR/XIONN(9,401),XIONV(9,401),XMASS(9)
      COMMON/MSS2/MASS(4),A(4)
      DATA XMASS/23.4164D-24,26.7616D-24,6.6904D-24,6*0.0/
      !... atomic masses and numbers for He+ and H+
      DATA MASS/6.6904E-24,1.6726E-24,26.7616E-24,0.0/,A/4,1,16,0/
      END
